Dynamics of tachyonic preheating after hybrid inflation 
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We study the instability of a scalar field at the end of hybrid inflation, using both analytical 
techniques and numerical simulations. We improve previous studies by taking the inflaton field fully 
into account, and show that the range of unstable modes depends sensitively on the velocity of the 
inflaton field, and thereby on the Hubble rate, at the end of inflation. If topological defects are 
formed, their number density is determined by the shortest unstable wavelength. Finally, we show 
that the oscillations of the inflaton field amplify the inhomogeneities in the energy density, leading 
to local symmetry restoration and faster thermalization. We believe this explains why tachyonic 
preheating is so effective in transferring energy away from the inflaton zero mode. 
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I. INTRODUCTION 

The dynamics of reheating after inflation have been 
studied intensively (l], |j ever since the first models of in- 
flation were formulated. In many cases, this process takes 
place non-perturbatively via a parametric resonance be- 
tween the inflaton and matter fields [[| |]. In this sce- 
nario, which is known as preheating, the long-wavelength 
matter modes are exponentially amplified to a very high 
effective temperature. This may have the consequence 
that certain symmetries, which are spontaneously bro- 
ken at the reheat temperature, arc temporarily restored 
during this non-thermal stage [|[ [f|. When the fields 
thermalize, the symmetry breaks again, leading possi- 
bly to formation of topological defects (?]]. In models 
of electroweak scale inflation, the period of non-thermal 
symmetry restoration may also explain the baryon asym- 
metry of the universe j|, |^, . 

However, it was argued more recently by Felder et 
al. [[HJ that in hybrid inflationary models, this paramet- 
ric resonance often fails to take place. In these models, 
inflation ends because a scalar field becomes unstable and 
its long-wavelength modes grow while the inflaton field 
rolls down the potential. The dynamics of this "tachy- 
onic" or spinodal instability have been studied by many 
authors y, g [If, [[§ |l| [l| [l§, using a range of 
different methods. 

Naively, the inflaton field would be expected to start 
oscillating once it reaches the minimum of the poten- 
tial. What Felder et al. [jllj observed was that the 
oscillations die away almost immediately, and coined 
the term "tachyonic preheating" to describe this phe- 
nomenon. This damping rules out the possibility of ordi- 
nary preheating, but the tachyonic instability itself still 
can lead to many similar effects, such as generation of 
baryon asymmetry [fill, [T7|, defect formation Jul and 
particle production |17|, |lq|. Nevertheless, the physical 
origin of the damping has remained unclear until now. 

In this paper, we use both detailed analytical calcu- 



lations and non-perturbative numerical simulations to 
study the dynamics of tachyonic preheating after hy- 
brid inflation. Analytical results show that the up- 
per limit, fc*, of the range of momenta amplified dur- 
ing the instability has a simple dependence on the ve- 
locity ip of the inflaton field at the instability point, 
fe* 1=3 (2mg(p) 1 ^ 3 13, 53. This is essentially equivalent to 
calculations of the number of topological defects formed 
in a phase transition ^lj . 

We confirm this result in our numerical simulations, in 
which we include also the inflaton as a dynamical field, 
by measuring the power spectrum directly and also by 
measuring the number of topological defects formed in 
the transition. Furthermore, we study the later evolu- 
tion of the system, and show that because the poten- 
tial in hybrid models is extremely flat above the critical 
value, the oscillations of the inflaton field form roughly 
spherical "hot spots" , inside which the energy density is 
very high and symmetry restoration can take place lo- 
cally. Depending on the values of the couplings, these 
hot spots may appear already during the first oscillation 
of the inflaton field. We believe this explains the strong 
dissipation reported in Ref. jll], where the inflaton field 
was observed to settle down in its minimum after only 
few oscillations. 

We carry out most of our simulations in two dimen- 
sions, because that allows us to treat a wider range of 
length scales in a single simulation, but in order to check 
that our results are not specific of two dimensions, we 
repeat some of the simulations in three dimensions. The 
breaking of global symmetries leads to the existence of 
massless Goldstone bosons, which are instead not present 
in the case of the breaking of local gauge symmetries. We 
recall that there is no experimental evidence for the exis- 
tence of such scalars whose couplings to the particles of 
the Standard Model of Particle Physics are severely con- 
strained by the present data. Therefore we study also 
the case in which the broken symmetry is a local gauge 
invariance. 
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The structure of the paper is the following: We start by 
reviewing the basic features of hybrid inflationary mod- 
els in Section II. In Section III, we study the instability 



of the scalar field analytically both in the approximation 
that the mass term changes instantaneously and in the 
more realistic case of hybrid inflation. We present the de- 
tails of our two-dimensional lattice simulations, together 
with the results for the power spectrum in Section [j"v| , 
and the results for the defect density in Section |v|. In 



Section VI , we show that the oscillations of the inflaton 
field lead t o "ho t spots" with local symmetry restoration. 
In Section VII, we present results of three-dimensional 



simulations showing that our findings apply to that case 
as well. 



II. HYBRID INFLATION 

In models of hybrid inflation [^2| , the inflaton field tp is 
coupled to another scalar field %, which becomes unstable 
at a certain critical value of tp = tp c . This causes the slow 
roll conditions to break down, and inflation ends. The 
simplest realization of this idea has the potential 



V(<p, X) = \m 2 y + i 5 W + J A (x 2 



(1) 



where tp and x are re al scalars. During inflation, the 
inflaton has a large value, tp 3> <p c = m/g, where we 
have defined m = A 1 / 2 ?;. Therefore the effective mass 
term of the x field, m 2 (</?) = —m 2 +g 2 tp 2 , is positive and 
the symmetry x ~~ * ~X is restored. Eventually, tp reaches 
tp c , and m^tp) becomes negative, implying that the % 
field becomes unstable, or tachyonic. 

The COBE measurements require that at tp = <^cobe, 
when the fluctuations in the cosmic microwave back- 
ground left the horizon, 



^cobe) 3/2 
M p 3 F'(<pcobe) 



5.3 x 10" 



(2) 



where M„ 



,,, - (8TTG)- 1 / 2 w 2.4 x 10 18 GeV is the reduced 
Planck mass. Assuming that the value of the potential 
and its derivative do not change significantly between 
'PCOBE and <p e , we find 



V'{tp 



COBE J 



■mltpc 



and ^(c^cobe) ~ Vq = m 4 /4A. Consequently, 
o 1 q to 5 



5.3 x 10- 4 8A 3 /2M3' 



(3) 



(4) 



which is typically much less than to 2 , and therefore we 
neglect it in the study of the dynamics of the system. If 
the slow roll condition V'(tp)Mp <C V(tp) is satisfied, we 
have 



XL 
3H 



V'Mp 



1 



Vo 



5.3 x 10- 4 V3Af 2 



(60U) 2 , (5) 



where H is the Hubble rate at the end of inflation. Even 
in GUT scale inflation, tp is relatively low, and in mod- 
els of electroweak scale inflation, we find tp ~ (10 _3 eV) 2 . 
Note, however, that if the potential has a more compli- 
cated shape, tp can be higher, such as in inverted hybrid 
models fl4|| . Therefore, we treat tp as a free parameter. 

The dynamics of the fields are described by the equa- 
tions of motion 

d 2 tp(t,x) = -2>Htp{t,x)-g\ 2 {t,x)tp{t,x), (6a) 
d 2 X (t,x) = -3Hx(t,S)+{m 2 -g 2 tp 2 (t,x)) X (t,S) 

-A x 3 (t,f). (6b) 

Let us first discuss the dynamics of the system in the 
simplest approximation, in which the field fluctuations 
are neglected altogether. The initial conditions are given 
by the field values at the end of inflation, when the infla- 
ton tp has just reached its critical value tp c and is rolling 
down the potential with velocity tp and the x field is still 
at rest in its minimum, 

tp(0, x) — tp c and tp(0,x) = —tp, 
X(0,x) = and x(0,f) = 0. (7) 

When tp rolls further down, the defect field x is dis- 
placed out of the symmetric phase. As far as it reacts 
much faster than the inflaton field, we can assume that 
the field x) reaches immediately the minimum at the 
given value of tp(t, x) and is given by 



X 2 (t,x) 



to 2 — g 2 tp 2 {t) 
A ' 



(8) 



Substituting %(t, x) in the potential V(x, f) we obtain 
the effective potential for the inflaton field: 

Vcfi (x^),p)^^m 2 tp 2 ~^. (9) 

At late times, tp oscillates around zero, and when the 
amplitude is small enough, we can write down an ap- 
proximate equation of motion 



^) + 3i^(i)~-w 2 ^), 



(10) 



where cj 2 = m 2 g 2 /X. Defining 7 = 3i//2, we can write 
the solution as 



tp(t)=Ae-^ cos (t^ 2 -^ 2 ). (11) 

Therefore the inflaton field shows a damped oscillatory 
behaviour with respect to t. 



III. INSTABILITY 

We are interested in understanding the dynamics of 
the defect field just after the inflaton field has crossed 
the critical point tp c . We take into account the fluctua- 
tions of the defect field. As the mass for the defect field 
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is negative, the modes at low momenta experience a fast 
growth. We ignore the fluctuations in the inflaton field, 
because at early times, as we shall mainly be interested 
in the case in which j 2 < A, the effects of the interac- 
tions between the fluctuations of the inflaton field and of 
the defect field can be neglected and, even if g 2 ~ A, the 
homogeneous mode of the tp field gives the dominant con- 
tribution. At late times, the linear approximation breaks 
down in any case and the results given in Sections 111A 
and IIIB are not valid anymore. 



A. Instantaneous quench 

Let us first discuss the instability in the case when the 
mass parameter of the \ field changes instantaneously 
from zero to —m 2 . This case has been studied before by 
many authors (see, for instance, Refs. jllL O, |ll| [TtJ). 
For simplicity, we ignore the expansion of the universe. 

In the linear approximation, the equation of motion for 
a Fourier mode x(t, k) of wave number k is 



as 



d 2 x(t,k) = (-k 2 + m 2 ) X (t,k). 



(12) 



Modes with |fc| < m become unstable and grow exponen- 
tially. 

In a quantum field theory, Eq. ( |l2| ) is valid as an op- 
erator equation, and in a classical field theory as a field 
equation. In the quantum theory, the initial state is the 
vacuum, and at tree level, it is completely described by 
the two-point functions of the fields, 



X*(k)x(k')) = ^-(2^ 3 5 3 (k-k% 
2\k\ 

7r*(fc)7r(fc')) = lM(2 7r ) 3 < y 3 (fc-fc'), (13) 



where tt = doX- It is important to note that as long as 
the equation of motion is linear, the operator nature of x 
only appears in these initial conditions, and therefore the 
classical field theory reproduces exactly the same time 
evolution, if the initial conditions are given by a Gaus- 
sian random field with the two-point function in Eq. (13). 
Thus, although we shall treat the theory as a classical sys- 
tem, the results of the linear approximation are identical 
to the quantum theory. 

Solving Eq. (|l2|), we find that the Fourier mode x(t, k) 
is given by 



x(t,k) = %(0, k) cosh ( ty m 2 — k 2 
tt(0, k) 



sinh [ty m 2 - k 2 . (14) 



V m 2 — k 2 

Therefore, the variance of the fluctuations (x 2 )(t) grows 



(X 2 )(t) = I fifX (x*(t,k) X (t,k')) 



(2tt) 3 (2tt) 3 



cxp(2mi), 



32?r 2 t 



sinh 2 (ty/ m 2 — A: 2 ) 
(15) 



where (x 2 )(0) is a divergent "vacuum" contribution. 

Another quantity we shall consider is the power spec- 
trum P(k), by which we mean the kinetic energy density 
of a given Fourier mode of the field Xi 

P(k){2Tr) 3 S 3 (k - k') = (<K*(t, k)w(t, k')). (16) 

Using Eqs. (O) and (fil), we find 



P{k) = tt + 777- sinh 2 tVm 2 - k 2 , (17) 
2 2k 

for k < m. This shows that modes with k<m are ampli- 
fied exponentially, or in other words, tachyonic preheat- 
ing produces a power spectrum with an effective cutoff 
scale fc* « m [[Lll Efj. 



B. Hybrid model 

The assumption that the mass parameter changes in- 
stantaneously makes analytical calculations easy, but it 
is not a very good approximation for the instability at 
the end of hybrid inflation. During inflation, the inflaton 
field ip is rolling down the potential slowly, and because 
it is only weakly coupled to the x field, the effective mass 
parameter m 2 e g(<p) changes typically very slowly. This 



5| using the linear ap- 
Of, Pll in the context of 



case has been studied in Ref. 
proximation, and also in Refs. 
defect formation. 

We assume that the fluctuations of the inflaton field ip 
are negligible, which means that we can describe it by its 
homogeneous part ip(t). Near the instability point, we 
can approximate 



(p(t) PH<p c - Ipt. 



(18) 



We are interested in the time-dependence of the fluc- 
tuations x(t, k) of the defect field, and at early times we 
can linearize Eq. (|6b|), 

d 2 oX (t,k) = {m 2 ~g 2 V 2 (t)-k 2 )x{t,k) 

~ (2m gi pt~k 2 ) X (t,k). (19) 

The solution for %(t, k) is given in terms of Airy functions: 

X(t, k) = c A (k)Ai fujt - K^j + c B {k)Bi (ut - , 

' (20) 
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where we have defined uj = (2mgip) 1 ^ 3 . 

This shows that each mode starts growing at the time 
t = tk = k 2 /uj 3 , when the effective mass becomes nega- 
tive m 2 , cS = —cu 3 t + k 2 < 0. We can express the coeffi- 
cients CA(k) and c B (k) in terms of the values of x(fc) and 
7r(fc) at the time t k as 



c A (k) 
c B (k) 



1 X(tk,k) n(t k ,k) 

2 Ai(0) wAi'(0) 

X(tk, k) _ ir(t k ,k) 
2V3 MO) uAi'(0) 



(21) 



where Ai(0) = 3- 2 / 3 r(2/3) _1 » 0.355 and Ai'(0) = 
-3- 1 / 3 r(l/3)- 1 « -0.259. 

In the very early time regime, each Fourier mode be- 
haves as 

X (t, k) = x(0, k) + (t- t fc )7r(0, k) + O ((t - t fc ) 3 ) , (22) 
which implies 



(x 2 )W-(x 2 )(o) 



o 

i 2 (u; 3 t) ; 



(D + l)(D + 3)(i? + 5)' 



(23) 



where D is the dimensionality of space and CVj 
2(47r)-' D / 2 r(£)/2)- 1 . In the physical case D = 3, 



<x 2 X*)-<x 2 >(o) 



6+4 



96tt 



2 ' 



(24) 



The spinodal value of the field x, on the other hand, 
grows linearly, 



Ddal 



(*) 



3A 



uj 3 t 
3A ' 



and therefore at early times, 

(x 2 )W-(x 2 )(0)«x s 2 pin odaiW- 



(25) 



(26) 



This result shows that no matter how small ip is, there 
is always a period after the instability, during which the 
linear approximation is valid. Eqs. (|24| ) and (26) seem to 
imply that the back reaction sets in when 



.3D 
spin 



32tt' 



1/3 



(27) 



but because Eq. (|24|) is a power series in ut, the expansion 
has broken down by that time, and instead, we have to 
consider the asymptotic behaviour at late times. 

At late times, Ai(uit — k 2 /uj 2 ) decreases exponen- 
tially while Bi(ujt — k 2 /uj 2 ) shows an exponential growth. 
Therefore, we neglect the first term in the right hand side 



of Eq. (p0[). The asymptotic expansion of Bi(z) for large 
values of z is given by 



Bi(z)~ -U-^expf-z 3 / 2 

\/7T V 3 



(28) 



The variance of the defect field x(t, x) at late times can 
therefore be computed as 



(x 2 )(t) 



" H d D k d D k' 



x*(t,k) x (t,k') 



(2tt)D (2k)D 

C D J dkk - 1 (c B {k)*c B {k) 
xBi(uj(t-t k )) 2 , 



(29) 



where we have ignored the constant vacuum contribution 
(X 2 )(0). Using Eqs. (|lf) and @, we have 



c B {ky CB {k)) = - 



1 



1 / k 2 
— [ 2.645 + 4.976— 

UJ 1 



(30) 



Therefore, we find 



<x 2 )W-^f 

Jo 



dkk 



D-2 



2.645 



y/u(t - t k ) 



y/u(t-t k ) 



(31) 



which we can expand around small k to obtain the 
asymptotic behaviour 



( X 2 )(t) ~ 2.645 I £r 
' uj 3 t 



C n „{D-l 



16tt V 2 

D-l 
2 

exp 



(a;*) 4 x 



%tf 2 

3 V 7 



In the physical case D — 3, this simplifies to 



(x 2 )W 



2.645 uj 



64tt £ 



exp 



; (ujt) 



3/2 



(32) 



(33) 



Using Eq. (^0|) , it is also straightforward to write down 
the expression for the power spectrum P[k), 



P{k) = 



u 2 [Ai'{uj{t - t k )) + ^Bi'(uj(t - t k )) 
8k Ai(0) 2 

k [Ai'{uj{t-t k ))-^ % Bi'{uj{t-t k )) 



8 Ai>(0) 2 
At late times, t t k , this behaves as 
fuH ( uj 2 k 2 \ 



p(k) 



24irk \0.355 2 0.259 2 J 



fM*-**))" 2 , (35) 
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but the high-A; modes for which t < tk, are still in vac- 
uum. Therefore, we conclude that that power spectrum 
has a n eff ective cutoff: The energy density in modes with 
k < V luH is exponentially high, but modes with k > VtoH 
are still in vacuum. 

The growth of the long- wavelength modes stops around 
the spinodal time t sp - m , which we can estimate by com- 
paring Eq. (||^) with Eq. (|25| ) and ignoring all factors of 
order one, 



1 



spin 



3 to 

7 ln — r 

oj \ 4 A 



3-D \ 2/3 



(36) 



and analogous equations apply to tp and r. We have used 
here the lattice Laplacian 

A/(f } = ~h ? + %) ~ 2f[s) + /(f " i)] ' (40) 



where i is a vector of length Sx in the direction i. 

The lattice approach can easily be used in an expand- 
ing universe, as well, by using conformal coordinates 
dr\ = dt/a and rescaled fields in which the expansion of 
the universe only appears as time-dependent mass terms 
m 2 , m 2 in the potential in Eq. ([!]), 



which agrees with Eq. (30) of Ref. ]iq| . In contrast to 
Eq. (p7|), the dependence on A is only logarithmic. In 
most cases, the logarithm can be ignored, and we find 
the simple result t sp - ln m lu^ 1 . 

Assuming that the back reaction does not change the 
overall shape of the power spectrum P(k) significantly, 
we can see that the final value of the cutoff k is 



3 uj 
— In 

4 A 



3-D 



1/3 



(37) 



Apart from the logarithmic factor, which is usually neg- 
ligible, this result agrees with Ref. [|l3| 

Ignoring the logarithmic factor and using the slow-roll 
condition (||), we can write Eq. (p7|) as 



V3 5 



1/3 



5.3 x 10- 4 8A 



M 2 



(38) 



Consequently, the instability becomes weaker if the en- 
ergy scale of inflation is low, and only non-relativistic 
particles are produced. 



IV. SIMULATIONS 

In order to study the dynamics of the instability nu- 
merically, we approximate the system with a classical 
field theory, which is believed to be a good approxima- 
tion when occupation numbers are large ]23[ |. At early 
times, the quantum correlations may play an important 
role (see, e.g., ^5|), and while it would be interesting 
to study their consequences, that is beyond the scope of 
this paper. 

For the scalar theory in Eq. (Q), the discretization of 
space and time is straightforward: We replace the deriva- 
tives by finite differences. It is convenient to define the 
fields x and if at the lattice sites and their time deriva- 
tives 7r and t at half way between two time steps. The 
discretized equations of motion for x and ir are 

n(t + 5t/2,x)-Tr(t -5t/2,x) = ^ _ dV{ X ,v) 



St 

X{t + St, x) - x (t, x) 
St 



dx 

n(t + 5t/2,x), (39) 



2l \ 2 2 

'•—■ma 



m (rj) 



a 
a 



(41) 



where rj is the conformal time, m is the physical mass and 
the double prime indicates the second derivative with re- 
spect to the conformal time. We assume that the universe 
is radiation dominated, in which case a" = 0. 

In order to approximate the quantum time evolution 
with the classical equations of motion, we use initial con- 
ditions that mimic the properties of the quantum vac- 
uum (^3). We do this by generating a Gaussian random 
field with the two-point functions in Eq. (|l3|). This re- 
produces the quantum tadpole diagram correctly and can 
be interpreted as a leading-order perturbative quantum 
correction. However, it has the drawback that it gener- 
ates radiative corrections. This is a major problem for 
the inflaton mass term m 2 ,, because it was supposed to 
be extremely small, but the radiative correction gives a 
contribution of the order Sm^ — 0(g 2 /Sx) in two spatial 
dimensions and 5m 2 = 0(g 2 /Sx 2 ) in three, which is rel- 
atively large. Following Ref. JhJ, we solve this problem 
by "renormalizing" the mass term: We use a negative 
non-zero mass term —Sm^ in the equations of motion, 
and choose its value in such a way that the radiative cor- 
rection cancels. In two dimensions, we need to compute 



Sm 



9 



4ir 2 Sx J 



d 2 w 



' sm 2 u> T 



sin 2 u> T 



0.32|-. (42) 
Sx 



The corresponding three-dimensional result [ic[| for the 
potential in Eq. @ is (5m 2 w 0.452 5 2 /<5a; 2 . 

To study the instability and growth of the long- 
wavelength x modes after the symmetry breakdown, 
we used two-dimensional simulations, because we could 
reach larger system sizes and longer s imu lation times that 
way. As we shall discuss in Section VII, we also carried 



out three-dimensional simulations to make sure that, al- 
though the details of the back reaction are dependent on 
the dimensionality, the overall qualitative picture is not. 

In the two-dimensional simulations, we chose x to be 
a real one-component field x- The lattice size was L — 
2016 2 , and in units in which the mass parameter was 
m 2 = 1, the lattice spacing was Sx — 0.5. We used two 
set of couplings, g = 10~ 3 and g — 10~ 4 , with A = 200g 2 
in both cases. In order to isolate the dependence on (p, 
we used a Minkowskian space with H = 0. We recall 
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FIG. 1: Power spectra of the x field measured at various 
times after the transition for g = 10 -4 , ip = 1.0. The solid 
and dotted lines show the results of the numerical simulations 
and of the analytical approximation in Eq. ([m|), respectively. 
From bottom to top, the lines correspond to the instants when 



<X > = 1, 10, 100, 1000 and 10000. 
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FIG. 2: Power spectra of the \ field at two different times 
after the transition for g = 10 -4 and various different values 
of (p. The dashed and solid lines corresponds to the instants 
when (x 2 ) = 100 and 10000, respectively, and the four differ- 
ent curves correspond to values tp = 0.1, 1, 10 and 100 from 
bottom to top. The dotted vertical lines show the cutoff scales 
k, calculated in Eq. ( |37j ) (ignoring the logarithmic factor). 



that ip and H are related by Eq. (Eh but we treat <p as a 
free parameter. 

At the start of the simulation, the inflaton held was 
given an homogeneous initial value ipi — 1.01<y5 c and a 
certain initial velocity ip. As discussed above, the initial 
condition for the x held was a Gaussian random field that 



satisfied the "quantum" two-point function. Eq. (|13|). We 
followed the time evolution of the field configuration by 
solving numerically the equations of motion, Eq. (|3S|). 
The perturbative mass counterterm, Eq. (|42|), ensured 
that the inflaton ip did not experience any significant ac- 
celeration or deceleration before it reached ip c . 

We measured the power spectrum P{k) of the x field 



defined in Eq. (16) at the instants when (x*){t) was equal 
to some positive power of 10. In Fig. |lj, we compare 
the simulated power spectra to the analytical result in 
Eq. ( |34| ) for the parameter values g — 1CP 4 and tp = 1.0. 
We can see that at early times, the linear approxima- 
tion works extremely well, apart from the very lowest 
momenta. However, these momenta correspond to wave- 
lengths which are of the order of the system size, and 
therefore this discrepancy may well be due to finite-size 
effects. Another possibility is that the fluctuations, small 
as they are, give rise to enough friction to damp down the 
exponential growth of the modes with lowest k. Both of 
these effects are artifacts of our numerical technique, and 
would therefore be absent in the true quantum theory. 

At later times, when (x 2 ) ~ 1000, a discrepancy ap- 
pears at high momenta, due to the back reaction. Fur- 
thermore, when (x 2 ) ~ 10000, we can see that the linear 
approximation has broken down. 

In Fig. ^|, we show the power spectra measured at the 
instants when (x 2 ) ~ 100 and 10000 for different val- 
ues of tp. At the former time, the evolution is still lin- 
ear and well described by the analytical formula, given 
in Eq. (p4[), while at the latter time, the non-linearities 
have set it. Nevertheless, the power spectrum retains its 
overall shape, with very high values at low momenta and 
a sharp cutoff around the same momentum scale as in 
the linear approximation. In particular, the back reac- 
tion does not change the cutoff scale k*, and therefore 
the prediction k* ss (2mgp>) 1 ^ 3 of the linear theory in 
Eq. (B7y remains valid. 



V. DEFECT FORMATION 

Depending on the number of components in the scalar 
field x, the potential V(x, <p) is invariant under discrete 
or global symmetries, which get spontaneously broken 
when the inflaton tp crosses its critical value <p c , if x ob- 
tains a non-zero expectation value. In general, this leads 
to the formation of topological defects (£6| . These defects 
may have several important effects on the later evolution 
of the universe |^7j , and in many cases the fact that they 
are formed can by itself rule out the model. On the other 
hand, defect formation is also important, because it is a 
convenient way of probing the dynamics of the transition. 

In a scalar field theory, defect formation can be un- 
derstood in terms of the Kibble mechanism (For 
a recent review, see Ref. f28| .) When the symmetry 
breaks, all directions on the vacuum manifold are equally 
probable, and whichever the system happens to choose 
is purely a matter of chance. However the correlation 
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length of the \ ne ld cannot be infinite at the time of the 
transition, and therefore the system cannot choose the 
same direction everywhere in space but only inside do- 
mains whose size is determined by the correlation length 
£ at the time of the transition. 

This argument was made more accurate by Zurek [ fl9| , 
who related £ to the rate of the transition and the criti- 
cal indices of the theory. The present case in which the 
defects are formed at the end of inflation is particularly 
simple, because, as the inflation has diluted away all the 
energy density, we can assume that the universe is at 
zero temperature. Therefore the equilibrium correlation 
length is given by the inverse of the mass, and near ip c 
we have simply 




1000 



2 2 

9 <P 



-2gmtpt 



-wH. 



(43) 



cty/dt 



Furthermore, the Lorentz invariance guarantees that the 
relaxation time of the field x IS equal to the correlation 



length, 



When ip approaches ip c , the equilibrium values of both 
£ and r diverge, and eventually x must fall out of equi- 
librium. This happens approximately at the time when 
the relaxation time t is equal to \t\. Thus, we have 



- 1 /2 =CJ -1. 



(44) 



This determines the typical distance between the topo- 
logical defects. 

Comparing Eq. Q) with Eq. @, we see that | k" 1 , 
and this is obviously no accident. Indeed, a power 
spectrum with a cutoff fc* means that the field changes 
strongly on distances larger than but is smooth at 

shorter distances. It is therefore clear that the correlation 
length is given by 1/A;*. Thus, we are simply using differ- 
ent language to describe the same physical phenomenon. 

We tested the prediction in Eq. ( fi"4| ) in the two- 
dimensional simulations described in Section IV . Because 



X is a real one-component field, the broken symmetry is 
Z 2 and the topological defects are domain walls. With 
the parameters we were using, Eq. ([II]) predicts that the 
typical distance between walls is £ » 0.79(gp)^ 1 ^ 3 . We 
can also estimate the total length of domain walls imme- 
diately after the phase transition to be 



FIG. 3: The ratio of the measured wall length Lwaii to the 
prediction in Eq. (^) as a function of the velocity (p. The solid 
lines are for g = 10 -4 and, from top to bottom, correspond 
to the times when (\ 2 ) = 10 3 , 10 4 , 10 5 and when (ip) = 0. 



The dashed line with open circles corresponds to g 
and (x 2 > = 10 3 . 



10" 




j prod 
^ wall 



L 2 /£« 1.28 x 10 6 ( 5 <^) 1/3 . 



(45) 



In our simulations, we measured the total length of 
domain walls L wa ii in the system by counting all pairs 
of neighbouring lattice sites where x has different signs 
and multiplying the result by the lattice spacing Sx. We 
carried out these measurements at various times after 
the transition, and the results are shown in Fig. || for 
different values of (p. In the figure we plot the ratio of the 
measured value to the predicted one given by Eq. (45), 
and therefore a horizontal line corresponds to a power 
law with the predicted exponent 1/3. Each data point is 
an average of eight different runs with different random 



FIG. 4: An example of a network of domain walls after the 
transition. Blue (dark grey) and red (medium grey) regions 
correspond to positive and negative values of x> an d x van- 
ishes at the green (light grey) domain walls. 



initial conditions, but we used the same set of eight initial 
conditions for all values of tp. 

Testing the prediction in Eq. (Ejj) is not entirely 
straightforward, because, as the domain wall network 
evolves with time, the measured value of £ wa ii depends 
on when the measurement was carried out. This choice 
is more or less arbitrary, and we have used two different 
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ways of determining the time of the measurement: Either 
when (x 2 ) crosses a certain threshold value, or when (tp) 
crosses zero. These choices can be interpreted as different 
ways of giving the quantity L wa u a precise definition. 

When (x 2 ) is small, domain walls are ill-defined ob- 
jects, because the sign of x can vary simply due to 
random fluctuations. Therefore the measurements at 
(x 2 ) = 10 3 give high values and do not agree with the 
predicted power law. Later, when (x 2 ) = fO 4 or 10 5 , 
the correspondence between the measurements and the 
prediction ([l5]) is good at high values of tp, apart from 
a constant factor of around five. Because the prediction 
was intended to give the order of magnitude only, we can 
conclude that the agreement is good. 

However, we can see that at very low values of tp, the 
predicted power law breaks down. The most likely expla- 
nation for this is that these extra domain walls are gen- 
erated by vacuum fluctuations, which are not strongly 
enough suppressed when (x 2 ) is relatively small. This 
means that the corresponding definition of L wa ii is not 
suitable for such slow transitions. Of course, we cannot 
rule out the possibility that the prediction ( f45|) breaks 
down at low tp in the full quantum theory, either. 

We also tested how L wall depends on the strength of 
the couplings. The open circles in Fig. || show L waii /L^^ 
in simulations with g — 10~ 3 measured at the time when 
(x 2 ) = 10 3 . This corresponds to (x 2 ) = 10 5 in the simu- 
lations with g = 10 , and we see that they agree reason- 
ably well, and therefore the dependence on g predicted 
by Eq. ( |45|) seems to be correct. Furthermore, it seems 
that the velocity tp below which the prediction breaks 
down is independent of g. 

Finally, we also measured L wa n at the time when tp 
reaches zero for the first time, which is typically much 
later than any of the other measurements. We can see 
that at high tp, the power-law prediction breaks down in 
this case. This all illustrates that measuring L wa n accu- 
rately is rather difficult, because the result depends on 
the time when the measurement is carried out. Never- 



theless, our results show that the prediction (45) works 
reasonably well in a wide range of velocities tp. 



VI. LOCAL SYMMETRY RESTORATION 

In ordinary thermal phase transitions, the defect net- 
work formed in the transition evolves mostly dissipa- 
tively, and the domain wall loops would simply collapse 
and disappear with time. However, in the case of hybrid 
inflation, the dynamics of the inflaton field may have sig- 
nificant effects on the later evolution of the system. 

Eq. ( pd| ) shows that if the x ne ld were completely ho- 
mogeneous, the inflaton would oscillate with roughly the 
frequency uj 2 ~ V^(0) = (g 2 /X)m 2 . However, the pres- 
ence of a defect network means that x is f ar from being 
homogeneous: It vanishes in the cores of the defects, and 
furthermore, when the defects annihilate, they release 
energy, which heats up the system locally and leads to 



inhomogeneities in x- I n this case, the frequency of tp 
becomes dependent on the position in space. 

In regions where the effective lu v is lower than the av- 
erage value, the tp field will start to lag behind, and a 
gradient energy starts to build up. We can study this 
effect in more detail in a simple model of one scalar field 
with an inhomogeneous mass term, with the equation of 
motion 

d 2 tp~ V 2 tp + M 2 (x)tp = 0, (46) 
where we adopt a Gaussian shape of width A for M 2 , 

l-eexp(-^j) 



M z {x) = Mq 



1 



(47) 



We assume that the spatial extent of the perturbation is 
large compared with the mass, A 3> M _1 . 

We assume that initially, the field has the value 
p(0, x) = tp and is at rest. At the spatial infinity, 
|af| — > oo, the field oscillates with a constant amplitude 



tp(t, oo) = tp cos 



/ M t 



(48) 



Because A is large, M{x) is slowly varying, and it 
makes sense to treat the spatial gradient term in Eq. ( |4"6| ) 
as a small perturbation. Therefore, we write 

tp(t,x)=tp (t,x)+Stp(t,x), (49) 

where tp satisfies the ordinary differential equation 

d 2 tp (t, x) + M 2 (x)tp {t, x) = 0, (50) 

and has the solution 

tpo(t,x) = tpcos(M(x)t). (51) 

In the equation for dtp, we ignore its gradient term and 
get 

d 2 Stp + M 2 (x)Stp = V 2 tp 



= —tp t 



t(VM) 2 cos(Aft) + V 2 Msin(Aft) (52) 



This ordinary differential equation has the solution 



Stp(t,x) 



tp t 



2(\7 In M) 2 M 2 t 2 sin(Aft) 



12M 

-3 (v 2 lnAf ) (Mtcos(Mt) ~ sin(Mi)) (53) 



At the origin x — 0, the gradient VlnM vanishes and 
Eq. ([53| ) simplifies to 



5tp(t, 0) = p 



W 



8A 2 A/ 1 - e 



[sin(Af t) - M t cos(M t)} , 



(54) 



where D is the dimensionality of space. This behaviour 
is shown in Fig. [5] for the parameter values D = 1, tp = 



FIG. 5: The amplification of the oscillations by an inhomo- 
geneous mass term. The dotted line shows the unperturbed 
oscillation, and the solid and dashed lines show the full nu- 
merical solution and the analytical approximation in Eq. (5j) 
for the field at the origin with the perturbed mass term (47) 
and the parameter values D = 1, tp = 1, e = 0.5, A = 10 and 

Mi = i. 



1, e = 0.5, A = 10 and Mq = 1, together with the 
full numerical solution of Eq. ( |46| ) . We can see that the 
amplitude of the oscillation is growing with time. 

In the case of hybrid inflation, this means that the 
amplitude of the oscillations of the tp field can eventually 
exceed tp c at certain localized points. Above that value, 
the potential becomes flat, and \tp\ can shoot to very high 
values. 

A somewhat similar phenomenon has been discussed 
by McDonald p^j , who showed using a linear analysis 
that the self-interaction of the inflaton field may lead to 
the formation of "condensate lumps" . In our case, the 
interaction with the matter field \ P ia Y s a more impor- 
tant role, because the infiaton-matter coupling is stronger 
than the inflaton self-coupling and because the instabil- 
ity at the end of inflation generates large inhomogeneities 
in the x field. 

If the damping of the oscillations caused by the expan- 
sion of the universe or by interactions with other fields 
is very weak, this overshoot may happen already when tp 
reaches its turning point for the first time. A very rough 
way of estimating how likely this is consists in calculating 
the damping of the amplitude during the first oscillation, 



<p(Q) 



exp 



3ttH\ 

I ~ exp 



3m 

1m v 



(55) 



Thus, if m <C gM p , this localized overshoot is likely to 
happen at the first turning point. 

Even if there is no overshoot at that time, the am- 
plitude of oscillations grows gradually at points where 
the frequency is lower, and \tp\ may therefore exceed 
ip c at some later time. We can see this effect in our 



FIG. 6: A snapshot of the energy density in a two- 
dimensional simulation a relatively long time after the phase 
transition. The symmetry has been locally restored in certain 
circular regions by the oscillations of the inflaton field ip. 



two-dimensional simulations, even when we take into ac- 
count the expansion of the universe with the Hubble rate 
H = 0.0024, which corresponds to m/gM p w 0.1. First, 
the domain wall network is evolving, with cusps and loop 
collapses releasing shock waves of energy. After a few 
oscillations of tp, however, the energy density and the 
amplitude of the oscillations of tp in certain localized re- 
gions start to grow. Eventually, \<p\ exceeds <p c in these 
regions, and the Z2 symmetry is locally restored. Be- 
cause the potential of ip is extremely flat above tp c , the 
force that pulls ip back towards its minimum, effectively 
disappears at this point. 

A snapshot of the energy density at this stage is shown 
in Fig. 0. The bright, almost circular disks are regions in 
which the energy density is very high. In this particular 
case, all these "hot spots" are located at a domain wall, 
but this is not always the case. In some simulations, they 
appeared in places where shock waves form cusps or loop 
collapses hit each other. 

We stress that for the existence of these hot spots, it 
is crucial that in hybrid inflationary models, the inflaton 
potential is extremely flat above the critical value. In a 
model with a convex potential, they would only appear 
as regions with a slightly higher amplitude. Even with 
the flat potential, the gradient energy pulls ip towards the 
origin, and therefore it is slowly oscillating around zero 
with a high amplitude. Of course, there is nothing that 
would stabilize these hot spots, and therefore they even- 
tually die away having radiated away their energy. It is, 
nevertheless, interesting to speculate whether they could 
have significant cosmological consequences, because the 
local energy density inside them is much higher than 
the reheat temperature. In more complicated models in 
which the inflaton tp is charged under a continuous global 
symmetry, the hot spots may, in fact, become stable Q- 
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VII. THREE-DIMENSIONAL SIMULATIONS 

A. Global theory 

In order to make sure that our findings are not arti- 
facts of the two-dimensional theory, we also simulated a 
theory with three spatial dimensions. Because domain 
walls are ruled out in cosmology, we used a complex field 
X, which means that the topological defects in the theory 
are cosmic strings. The potential of the theory is 

Vfa X) = \™ 2 y + .9 2 |xl V + A (|x| 2 -v 2 ) 2 , (56) 

and the continuum equations of motion in conformal co- 
ordinates are |HJ 

9 V = didiv - 2g 2 \ X \ V (57a) 
d 2 Q X = d^X + (2A« V + g 2 tp 2 ) X - 2A| X | 2 X . (57b) 





(a) 



FIG. 7: Field configurations in the global run with couplings 
g = 0.01, A = 0.01 and m = 10 15 GeV. 

(a) The transparent surface shows the isosurface tp = —tp c at 
time t = 224, soon after the spatial average of tp reached its 
first turning point. The strings correspond to the isosurface 
\x\ 2 = f 2 /10 at the time when 

(b) The transparent surface shows the isosurface |x| 2 = v 2 /W 



for the first time. 

Ix! : 

at time t = 450, and the blue (dark) and yellow (light) opaque 
surfaces show the isosurfaces tp = tp c and tp — ~ip c at the same 
time. 



The time derivatives are with respect to the conformal 
time rj, and we have assumed that the universe is ra- 
diation dominated so that the scale factor behaves as 
a = 1 + Hrj. 

The discretization of the field equations was carried 
out in the same way as in Sect. [V, with lattice size 256 3 , 



lattice spacing Sx = 1.0 and time step St = 0.05 in units 
with m 2 = 2Xv 2 = 1. 

To compare with the results in Ref. Jll]], we carried 
out a simulation with the same coupling values g — 0.01, 
A = 0.01, m = 10 15 GeV. In our units, the Hubble rate 
becomes H — 0.0012. The initial velocity of the inflaton 
field was <p = —0.00814, and its initial value was such 
that it reaches tp c at time t = 10, i.e., <p(0) = 100.0814. 
As in Sect. |y|, the fields also had Gaussian "quantum" 
fluctuations in the initial configuration. 

In this case m/gM p as 0.04 -C 1, and indeed, the in- 
flaton tp overshoots — tp c in a significant fraction of the 
space, as shown in Fig. ^a. We have also plotted in 
the same figure an isosurface of |x| 2 at the time when 
tp crosses zero. It shows the cosmic strings formed in the 
transition. The regions with tp < ~tp c follow closely the 
shape of the cosmic string network, in agreement with the 
discussion in Sect. VI. In these regions, the inflaton field 
starts to lag behind and this effect destroys the homo- 
geneity of the oscillations. Indeed, Fig. || shows that the 
oscillations of the spatial average of tp die away during 
the first period, as observed earlier in Ref. 111. In fact 



, a 
at 



sign of this overshoot can be seen in Fig. 3 of Ref. JL1[ 
t = 130, the histogram has two peaks, one in the broken 
phase, corresponding to the bulk of the space, and one 
at tp < —tp c . 

Even though the spatial average does not oscillate, the 
space is full of oscillating hot spots even at late times, as 
shown in Fig. ^b. The energy is concentrated in these hot 
spots, and although tp oscillates with a high amplitude 




■</7 m=10'°GeV 

m=3x10 ,5 GeV 

m=3x10 l5 GeV, gauge theory 
m=6x10 ,5 GeV 
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time t 
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FIG. 8: Time evolution of the spatial average of the inflaton 
field tp. The origin of the time axis has been shifted in each 
curve so that the curves coincide when tp first rolls down the 
potential. The x's and +'s shows the times when the \ and 
tp configurations were plotted in Figs. (?], and [To| 



inside each hot spot, the phases are uncorrelated and do 
not contribute to the spatial average of tp. 

Eq. ([55]) predicts that if the mass scale m is increased, 
the probability of an overshoot decreases, and the damp- 
ing is therefore weaker. We studied this by repeating the 
simulation with m = 3 x 10 15 GeV, which corresponds 
to H = 0.00361. The initial values of the inflaton field 
and its velocity were tp(0) = 100.732 and tp = -0.0732. 
For these parameters, m/gM p ss 0.12, and Fig. ||a con- 
firms that tp exceeds —tp c only in some small, isolated 
regions. The spatial average of tp is shown Fig. |8|, and 
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FIG. 9: Field configurations in the global run with couplings 
g = 0.01, A = 0.01 and m = 3 x 10 15 GeV. 

(a) The transparent surface shows the isosurface <p = — <p c at 
time t = 130, soon after the spatial average of p reached its 
first turning point. The strings correspond to the isosurface 
|x| 2 = u 2 /10 at the time when <p = for the first time. 

(b) The transparent surface shows the isosurface |x| 2 = v 2 /W 
at time t = 270, and the blue (dark) and yellow (light) opaque 
surfaces show the isosurfaces <p — p c and ip = — tp c at the same 
time. 



during the first oscillation, the damping is largely due to 
the expansion of the universe [see Eq. (p7T|)] . 

Nevertheless, as discussed in Sect. [Vl|, inhomogeneities 
in the energy density amplify the oscillations locally. 
Eventually, hot spots are formed, and at late times, the 
system is again full of them, as shown in Fig. ||b. This 
leads to the strong damping at late times shown in Fig. [|. 

In Fig. ||, we also show the time evolution of the spa- 
tial average of ip in a run with m — 6 x 10 15 GeV, and 
in agreement with the arguments in Sect, [vj, the damp- 
ing is even slower in this case. Higher m corresponds to 
higher Hubble rate H, and therefore the fact that damp- 
ing becomes weaker when m increases is opposite to what 
the naive tree-level treatment in Eq. (O) predicts. 



B. Gauge theory 

Spontaneous breakdown of a continuous global symme- 
try gives rise to massless Goldstone bosons, which have 
not been observed in nature, and therefore it is interest- 
ing to consider the case in which the broken symmetry 
is a local gauge invariance. We need to couple the field 
X to an Abelian gauge field by replacing the derivatives 
in Eq. ( |57b| ) by covariant derivatives Di = di +ieAi, and 
introducing the equations of motion for the gauge field 



doEi 
d % Ei 



djF i:j + 2e Imx* Ax, 



2e Imx*<9 x, 



(58a) 
(58b) 



where Ei = —doAi, F^ = diAj 



The initial condi- 



tions for the gauge field are analogous to Eq. (|13|), e xcep t 
that they have to satisfy the Gauss constraint Eq. ( |58b[ ) . 

We carried out the discretization of the gauge field 
equations in the standard way as discussed, say, in 




FIG. 10: A field configuration in the gauge theory with 
couplings g = 0.01, A = 0.01 and m = 3 x 10 15 GeV. The 
transparent surface shows the isosurface |x| 2 = u 2 /10 at time 
t = 270, and the blue (dark) and yellow (light) opaque sur- 
faces show the isosurfaces p — ip c and tp — —p c at the same 
time. 



Ref. |3^] , with the scalar fields defined at the lattice sites 
and the gauge field on links between the sites. The scalar 
couplings were g = 0.01, A = 0.01 and m = 3x 10 15 GeV, 
and the gauge coupling was e = 0.1. For the Hubble rate, 
these parameters correspond to the value H = 0.00361. 
The initial values of the inflaton field and its velocity 
were <p = 100.732 and (p = -0.0732. 

From Fig. || we can see that the gauge field does not 
change the dynamics significantly. Therefore, the exis- 
tence of Goldstone modes is not crucial for the damping. 

In thermal phase transitions the gauge field plays a 
crucial role in defect formation 28 , 33) , but its effect de- 
pends on the temperature, and as the transition takes 
place at zero temperature, the gauge fields can be ne- 
glected. Indeed, the total length of string immediately 
after the transition is similar f or g lobal and local theo- 
ries. Therefore we expect Eq. (J44|) to be valid in gauge 
field theories, as well. Because global strings have long- 
range interactions while local strings do not, the string 
network decays much faster in the global case. 

The result that Eq. ( |44| ) applies to gauge theories, sup- 
ports the analytic estimates given in Ref. Jl3] for baryon 
density in a scenario in which the baryon asymmetry is 
generated by unwinding electroweak textures formed by 
the Kibble mechanism at the end of electroweak-scale 
hybrid inflation. The predicted number density of these 
"knots" is ph ~ £ -3 ~ ^* ~ 2mg(p. Once formed, they 
can decay either by changing the topology of the Higgs 
or the SU(2) gauge field configuration, and as pointed 
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out by Turok and Zadrozny |34], [55|, CP violation bi- 
ases the decay of the winding configurations. Thus, if 
CP violation is present, this mechanism leads to baryon 
asymmetry. For a reliable calculation of the baryon den- 
sity, one has to understand the dynamics of the gauge 
field in the presence of CP violation. 

The further evolution of the system confirms that the 
hot spots with local symmetry restoration appear also in 
the gauge theory. Fig. |n] corresponds to Figs. |?]b and 
Figs. ||b, and shows a very similar set of hot spots. 

VIII. CONCLUSIONS 

In this paper, we have studied tachyonic preheating 
in detail using both analytical and numerical techniques, 
and taking the full dynamics of the inflaton field into ac- 
count. We simulated the non-perturbative dynamics af- 
ter the spinodal time in two and three dimensional scalar 
field theories and a three-dimensional gauge field theory, 
and showed that at early times, the instability is well de- 
scribed by a linear approximation. When non-linearities 
set in, the long-wavelength fluctuations form topological 
defects, and our results show that their number density 
is related to the cutoff scale fc* in the same way as in 
thermal phase transitions. 

We also studied the dynamics of the system at later 



times, and showed that oscillations of the inflaton field 
coupled to the inhomogeneities of the matter field lead to 
hot spots of high energy density, inside of which the sym- 
metry can become temporarily restored. Depending on 
the coupling constants, this effect may be very strong, 
in which case it damps down the homogeneous oscilla- 
tions of the inflaton field. We believe this explains |3(| 
the high effectiveness of tachyonic preheating observed 
in Ref. ]ll[ . Further work is required for relating the 
density and sizes of these hot spots to the parameters 
of the model, and for understanding their other possible 
consequences. 
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